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Using  rer.ormalization-grou?  methods  and  the  postulated  equivalence  between  the  mental-range 
structures  o(  turbulent  flows  satisfying  initial  and  boundary  conditions  and  of  flows  driven  by  a  ran¬ 
dom  force,  we  evaluate  the  Kolmogorov  constant  (1.617)  and  Batchelor  constant  (1.161),  skewness 
factor  (0.4878).  power-law  exponent  (1.3307)  for  the  decay  of  homogeneous  turbulence,  turbulent 
Prandti  number  (0.7179),  and  von  Karman  constant  (0.372).  This  renormaliaation-group  tech¬ 
nique  has  also  been  used  to  derive  turbulent  transport  models. 

PACS  numbers:  47.25.-c 


The  direct  interaction  approximation  (DIA),  due  to 
Kraichnan.1  was  the  first  fieid-theoreticai  approach  to 
the  theory  of  turbulence.  Formulated  in  terms  of  the 
Dyson  equation,  the  DlA  is  characterized  as  the 
lowest-order  approximation  which  includes  nonlinear 
corrections  to  the  propagator  for  the  mode  v(k.w).  It 
was  shown1  that,  in  the  inertial  range,  the  DIA  gives 
the  energy  spectrum  £(£)«  k~in.  This  result  contra¬ 
dicts  both  experimental  data  and  the  Kolmogorov  the¬ 
ory  of  turbulence  which  gives  £(£)«  k~5,i,  perhaps 
with  small  corrections  due  to  intermittency. 

The  source  of  this  discrepancy  between  the  DIA  and 
the  Kolmogorov  theory  has  long  been  understood.2 
The  DlA  does  not  distinguish  between  dynamic  and 
kinematic  interactions  between  eddies  of  widely  sepa¬ 
rated  length  scales.  Small  eddies  are  convected  by 
large  eddies  in  a  purely  kinematic  way  which  should 
not  lead  to  energy  redistribution  between  scales.  The 
spurious  effect  of  large-scale  convection  on  small 
scales  has  been  removed  from  the  DIA  by  use  of  a 
Lagrangean  description  of  the  flow.  This  Lagrangean- 


history  direct  interaction  approximation3  (LHDIA) 
leads  to  the  Koiomogorov  5/3-energy  spectrum  with 
the  Kolmogorov  constant  Ck-1.77  [see  (11)  below] 
which  is  in  reasonable  agreement  with  experiment.4 
However,  application  of  the  LHDIA  to  the  problem  of 
turbulent  diffusion  of  a  passive  scalar  does  not  lead  to 
quantitative  agreement  with  experimental  data:  The 
turbulent  Prandti  number  P,  calculated4  from  the 
LHDIA  is  roughly  0.14,  much  smaller  than  the  experi¬ 
mentally  observed  P,  *»  0.7-0.9. 

In  1977  Forster.  Nelson,  and  Stephen5  used  dynamic 
renormali2ation-group  (RG)  methods,  originally  de¬ 
veloped  for  the  description  of  the  dynamics  of  critical 
phenomena,4  to  derive  velocity'  correlations  generated 
by  the  Navier-Stokes  equation  with  a  random-force 
term.  The  ideas  expressed  in  Ref.  5  have  been  used 
by  others  in  the  context  of  hydrodynamic  tur¬ 
bulence.7-10  The  problem  is  formulated  as  follows: 
Consider  the  (/-dimensional  space-time  Fourier- 
transformed  Navier-Stokes  equation  for  incompressi¬ 
ble  flow. 


vlU)^G0Mk)-{ik0G0Plim,^)JvM(q)vl,U-q)ddq/(2tr)^\  (1) 

where  the  zero-mean  Gaussian  random  force  f(k,a>)  is  determined  by  its  correlation  function 

(j;(k.cu)/J(<,Ui'))  ~(2T)^*,2D0k-'Pu(k)5(k^k').  (2) 

Here 


G°-  ( -  uit  +  vok2)'1,  PtJ( k)  -8y-  k,kjl k1. 


£yjk(k)  -  kkPj(k)  +  kjPj'ik), 


(3) 


k  —  (k.ai),  v0  is  the  kinematic  viscosity,  X0-l,  and 
the  constant  y  >  —  2.  The  problem  (1M3)  is  formu¬ 
lated  on  the  interval  0  <  k  .\g  and  -«  <  u  <  », 
where  Ao  is  a  wave  number  beyond  the  dissipation 
wave  number  at  which  substantial  modal  excitations 
cease.  The  parameter  D0,  which  determines  the  inten¬ 
sity  of  the  random  force,  is  discussed  below. 

The  RG  procedure  consists  of  the  elimination  of 
modes  v>(k)  with  wave  vectors  satisfying  A0e-' 


<  k  <  A  from  the  equations  of  motion  for  the  modes 
v<(k)  with  wave  vectors  from  the  interval 
0  <  k  <  A0e-r.  At  this  stage,  kinematic  interactions 
are  excluded  by  construction  and  one  can  expect  phys¬ 
ically  meaningful  results  in  the  limit  k—  0.  Detaiis  of 
this  RG  procedure  are  given  elsewhere.5,11 

The  RG  scale-eliminauon  procedure  gives  a  correc¬ 
tion  to  the  bare  viscosity  v0  in  terms  of  an  effective 
viscosity  which  takes  into  account  the  effect  of  the 
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eliminated  modes.  The  result  is 

vr  —  Vq[1  •*"  .-ltf.Vg(  d‘f—  1  )/  «  !, 

where  «  — 4—y -  d,  Ad-  AdSji'(l* K  and 

-  1  <r-d-t  -  C-)J/: 

J~  2  d(d- 2)  •  f(i-j)  ' 


(4) 

(5) 


The  dimensionless  expansion  parameter  \0  (which  is  a 
Reynolds  number)  is  defined  as  A5  —  Dyv\ Ag.  As  we 
shall  see  below,  the  choice  of  /  —  d  recovers  the  Kol¬ 
mogorov  scaling  in  the  inertial  range. 

By  variation  of  the  cutoff  A  ( /•)  -  Ag e ” '  we  derive 
differential-recursion  relations  for  \(r)  -  {Drjvir)* 
x  A  ( r) 1/2  and  v(r) 

2^2.  —  Adv(r)\}{r),  --7-  — A*(«  — 3  Ad\~).  (6) 

dr  dr 

The  solutions  to  (6)  are 

l ( r)  -  \0e4'/2 1  1  +  3^(  e"- 1 )/« ]  “ 

v(r)  —  v0[  1  «f  3Ad\o(e,r~  1  )/<  1 1/3 . 

In  the  limit  r—  «  the  coupling  parameter  \  (which  is 
an  effective  Reynolds  number)  goes  to  the  fixed  point 

A.-(«/3  Ad)^‘ 

and 

Eliminating  all  modes  with  q  >  k,  we  set  A  -  kind  ob¬ 
tain 

v(k)-(±A<2i)uik‘''3  (7) 

-0.4217(2D,V(2w)'l1/3lt’4/3 

when  y  —  d  —  3.  The  coefficient  Ad  is  computed  from 
(5)  in  the  lowest  order  of  «  expansion  («  —  0);  thus 
/i*  -  0.2  in  the  three-dimensional  case  J  —  3. 

The  energy  spectrum  can  be  calculated  to  lowest  or¬ 
der  in  «  from  the  equation  v(k)  -  G(k)f(k),  where 
the  propagator  G(k)  is  evaluated  with  the  ^-depen¬ 
dent  viscosity  (7).  The  result  is 

E(  k)  -  l.l86(2A>S*,(2»)rf!“3A:":/3.  (8) 

Thus  the  renormalization-group  procedure  applied  to 
randomly  stirred  fluid  gives  the  Kolmogorov  spectrum 
in  the  case  y  —  d. 

In  order  to  complete  the  analysis,  it  is  necessary  to 
relate  the  parameter  Di)  to  observables.  Consider  a 
fluid  described  by  the  Navier-Stoices  equation 

^■  +  (v'V)*«--7jj*»)71v,  (9) 

of  p 

subject  to  initial  and  boundary  conditions.  We  assume 
that  strongly  turbulent  fluid  is  characterized  in  the 


ir.eruai  range  of  scales  by  statistically  universal  scahr.a 
iaws  (Koimogorov  spectrum,  etc.)  wrnch  are  indepen¬ 
dent  of  initial  and  boundary  conditions.  Thus,  the  sys¬ 
tem  in  the  universal  regime  can  be  described  by  equa¬ 
tions  of  motion  which  do  not  involve  any  particular  in¬ 
itial  and  boundary  conditions:  (1)  and  '2),  for  exam- 
pie.  provided  that  the  random  force  in  (1)  ar.d  (2)  is 
chosen  in  such  a  way  that  it  generates  velocity  fluctua¬ 
tions  which  are  statistically  equivalent  to  the  solutions 
of  (9)  subject  to  initial  and  boundary  conditions.  In 
other  words,  to  describe  the  fluid  in  the  inertial  range 
we  may  replace  (9)  with  the  corresponding  system 
(1),(2)  with  a  properly  chosen  force.  In  this  case,  it 
has  been  shown12  that  if  we  assume  that  solutions  of 
Eq.  (9)  in  the  inertial  range  scale  as 

v(lt)-i VSUik~^  (10) 


and 

£U)-CK.S2/3ATJ/3.  (11) 

then  energy  balance  in  analytical  turbulence  theory  re¬ 
quires  that  the  Kolmogorov  constant  CK  in  (11)  and 
the  parameter  N\n  (10)  be  related  as 

NtCi  -0.1904. 

Here  5  is  the  rate  of  energy  dissipation  in  the  fluid. 
Demanding  the  equivalence  of  (7)  and  (8)  with  (10) 
and  (11)  in  the  inertial  range  gives 

1.594  5,  (12) 

so  that  Cic“  1.617. 

A  similar  RG  procedure11  applied  to  the  equation  of 
a  passive  scalar  gives  the  result  that  the  turbulent 
Prandtl  number  P,  in  the  case  y  -  d  -  3  is 


1.3929, 


so  that  ^,  —  0.7179.  The  Batchelor  constant  Cg*  is  de¬ 
fined  by  the  inertial-range  scalar  fluctuation  spectrum. 
Using  energy  balance  in  terms  of  the  /^-dependent  vis¬ 
cosity  at  the  fixed  point,  we  find11  Cg,  -  CKP,  so  that 
C9l-  1.161.  Another  calculation*3  of  Cg,,  based  on  an 
RG-modified  version  of  the  direct-interaction  approxi¬ 
mation.  gives  the  same  result.  The  results  for  the  tur¬ 
bulent  Prandtl  number  and  the  Batchelor  constant  are 
in  close  agreement  with  experimental  data.14 

The  renormalization-group  procedure  can  also  be 
used  for  deriving  averages  of  different  nonlinear 
operators  over  the  fluctuating  velocity  field.11  For  ex¬ 
ample,  the  skewness  factor,  wtuch  is  a  dimensionless 
measure  of  nonlinear  transfer,  is  defined  as 


(Ov./3x,)3)  _  a 

(Ov,/^*,)2)3'2  "  B>n' 


where  (...)  denotes  average  over  the  fluctuating  velocity  field,  and 

A  -  ((oui/a.v,)3)  -  -  if  Q\P\(k  -  q  -  p){vl(q)vl(  p)v.(k  -  q  -  p)  ddq  ddpi  ( 2~ 

in  the  limit  k  —  0.  Decomposing  the  velocity  field  into  the  components  v*  and  vK  and  eliminating  small  scales 
using  the  forced  Navier-Stokes  equation  (1).(2),  we  find,  in  the  lowest  order  in  the  t  expansion,  that11 

AK  -  -  i  f  <t\P](k-q-p])vl<(q)v{c  ( p)vf  (k  -  q  - p)  ddq  ddp/ (2ir)ld+2 

--*iot2^oV(2-)rf]5/^A2 

in  the  limit  k— *  0  (r—  oo).  The  same  procedure  appiied  to  evaluation  of  3 in  (13)  gives 
D<  1  2iW(2r)' 

B  "  20  v 

in  the  limit  k —  0  (r—  »).  Thus 


A<  ID^SJdsrY 


■0.4878 


when  calculated  at  the  fixed  point  of  the  RC  calcula¬ 
tion.  Since  $  <  ( r)  does  not  depend  on  r  in  the  limit 
r—  oo,  we  assume  that  (14)  holds  everywhere  in  the 
inertial  range,  and  so  #  -  0.4878.  It  should  also  be 
noted  that  the  same  RG  procedure  gives  the  exact 
result  4*  —  0  in  the  rwo-dimensionai  case  d  - 2. 

Another  important  relation  can  be  derived  from  the 
Kolmogorov  energy  spectrum  and  formula  (7)  for  the 
turbulent  viscosity.  It  can  be  checked  readily  that  the 
total  kinetic  energy  K  in  the  system  is  K - 1.1955/ 
v\2,  where  A  is  the  wave  vector  corresponding  to  the 
integral  scaie  of  turbulence.  Combining  this  reiadon 
with  (7)  and  (12)  we  derive  a  relation  between  v, 
kinetic  energy  K,  and  the  mean  dissipation  rate  6 , 
namely,  v  —  0.0837AT2/  S . 

The  RG  procedure  can  be  used  to  evaluate  each 
term  of  the  equations  of  motion  for  kinetic  energy  and 
dissipation  rate.  This  leads  to  a  so  called  K  -  S  model 
of  turbulence.  It  can  be  shown11  that  this  RG  model 
implies  that  isotropic  turbulence  decays  as 
K<*  (/  —  to)"lJ307  which  is  close  to  the  experimental 
data1*  and  recent  results  of  direct  numerical  simula¬ 
tions.13  The  same  model,  which  does  not  involve  any 
experimentally  adjustable  parameters,  gives  the  von 
Karman  constant11  * —  0.372  for  the  logarithmic  velo¬ 
city  profile. 

The  good  agreement  of  the  RG-predicted  constants 
(C*.  Cb«»  S ,  x)  wuh  experimental  data  is  to  some 
extent  surprising  since  the  RG  procedure  does  not 
take  into  account  local  interactions  between  eddies  of 
similar  size.  However,  it  has  been  pointed  out9  that 
the  ratio  of  time  constants  which  correspond  to  nonlo¬ 
cal  and  local  interactions  is  0(«1/:).  Thus,  local  in¬ 
teractions  are  weak  if  <  is  assumed  small.  It  remains  to 
be  explained  why  the  lowesi-order  truncation  of  the 


RG  expansion  in  powers  of  €“4  works  so  well. 
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Abstract  -The  expression  for  turbulent  Prandtl  number  obtained  from  the  renormalization  croup 
procedure  is  used  to  describe  the  process  of  heat  transfer  in  turbulent  pipe  floss  The  results  are  in  a  good 
agreement  with  experimental  data  oxer  the  entire  range  of  experimentally  accessible  Prandtl  numbers. 

10" 2  <  <r„  <  I0\ 


1.  INTRODUCTION 


The  problem  of  heat  conduction  in  turbulent  flows 
has  been  under  intensive  study  for  more  than  half  a 
century.  Experimental  data  on  velocity  and 
temperature  distributions  have  suggested  many  semi- 
empirical  theories  to  describe  the  basic  properties  of 
the  phenomenon. 

It  has  long  been  realised  that,  if  the  Reynolds 
number  is  large  enough  and  the  Prandtl  number 
a0  I’o/ix'o  is  not  too  small,  the  molecular  difTusivity 
k0  does  not  play  any  role  in  the  process  of  heat 
conduction  or  diffusion  in  turbulence.  In  this  case,  the 
temperature  and  velocity  distributions  have  similar 
behavior  in  the  wall  region,  both  obeying  the 
logarithmic  law  with  the  temperature  profile 


and  semi-empirical  relations  to  describe  turbulent 
heat  transfer  across  a  wide  range  of  Prandtl  and 
Reynolds  numbers.  More  than  30  formulae  of  this 
kind  have  been  reviewed  by  Reynolds  [1]  in  1975.  In 
1979.  Gori  ef  al.  [2]  concluded  that  there  is  no  general 
way  to  describe  turbulent  heat  transfer  in  low-Prandtl- 
number  fluids  for  a  wide  range  of  Re.  They  suggested 
the  following  formula  for  the  turbulent  Prandtl 


number  o,urtl  when  Re  <  1.7  x  10s : 


O.OURe0-4*^* 


x  ;i-exp[-<0.014Re0*5ffg:r'];  (Ii 


as  proposed  by  Aoki  [3]  or 


( 1  +  IOOPe"05)^  1  +  120Re" 


-0.15] 


’fcpfpU, 


(In  y  +  C). 


as  proposed  by  Reynolds  [1],  The  formula 


Here  <  )  denotes  a  horizontal  average,  y  is  the 
distance  to  the  wall,  q  denotes  the  constant  heat  flux 
and  cp  and  u,  are  the  heat  capacity  and  friction 
velocity,  respectively.  The  Von  Karman  constant 


proposed  by  Jischa  and  Rieke  [4],  was  suggested  to 
represent  the  Reynolds  number  range  1.7  x  105  <  Re 


<  2.6 x  105;  the  constant  a,., 


0  85  was  used  for 


b  =  0.4  and  <r,urt)  =  vlllrt,/Kwrt,  is  the  ratio  of  turbulent 


viscosity  to  turbulent  heat  conductivity.  According  to 
the  well-known  Prandtl-Reynolds-Colbum  analogy, 
the  turbulent  Prandtl  number  is  nearly  a  universal 
constant:  erturt>  =*  0. 7-0.9. 

In  the  limiting  case  of  small  Prandtl  number,  the 
molecular  difTusivity  k0  cannot  be  neglected  and  the 
simple  analogy  between  temperature  and  velocity 
distributions  does  not  work.  It  is  clear,  however,  that 
as  c o  -•  0.  the  Nusselt  number  Nu  [defined  below  as 
the  dimensionless  (based  on  the  bulk  temperature,  see 
equation  (28))  heat  flux]  satisfies  Nu  %  const.  It  is 
known  from  experiments  that  Nu  %  6. 8-7.0  in  flows 
with  constant  heat  flux  through  the  wall  while  Nu  is 
somewhat  smaller  in  flows  with  constant  wall 
temperature.  To  the  best  of  our  knowledge,  there  is  no 
satisfactory  theory  describing  heat  conductivity  in 
turbulent  flow  with  low  Prandtl  number. 

Many  attempts  have  been  made  to  find  empirical 


Re  >  2.6x  105.  When  relations  (1H3)  are  used  to 
predict  the  mean  temperature  field,  they  give 
reasonably  accurate  predictionsof  the  Nusselt  number 
(which  is  related  to  the  wall  gradient  of  the 
temperature  profile).  However,  the  full  temperature 
profiles  predicted  on  the  basis  of  expressions  i  Ih3> 
were  less  satisfactory. 

In  this  work  we  apply  a  formula  for  the  turbulent 
Prandtl  number  derived  by  Yakhot  and  Orszag  [5]  to 
describe  heat  transfer  in  pipe  flows.  It  will  be  shown  in 
Section  3  that  the  proposed  relation  between  turbulent 
viscosity  and  turbulent  heat  conductivity  gives 
accurate  predictions  of  both  Nusselt  number  ami 
temperature  distributions  across  an  extremely  wide 
range  of  Prandtl  and  Nusselt  numbers. 


2.  FORMULAE  FOR  TURBULENT 
PRANDTL  NUMBER 


♦  Department  of  Mechanical  Engineering.  University  of 
Ben  Gurion  of  the  Negev,  Beer-Sheva  84105.  Israel 


Here  we  present  some  of  the  basic  ideas  leading  to 
an  expression  for  the  turbulent  Prandtl  number  The 
main  steps  of  the  renormalization  group  procedure  are 
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NOMENCLATURE 


A  Van  Driest  parameter 

B  proportionality  coefficient  in 

.r-coordinate  dependence  of 
temperature,  equation  (20) 

At  geometric  factor,  equation  (11) 

C  constant  in  the  temperature  profile, 

equation  (30) 

Ch  heat  transfer  coefficient 

D  pipe  diameter 

G°  bare  propagator  for  velocity 

Nu  Nusselt  number  (based  on  the  bulk 

temperature),  C^a0Re 
Pe  Peclet  number,  a0Re 

R  radius  of  pipe 

Re  Reynolds  number  based  on  the  pipe 
diameter,  u„D/v0 

R ,  Reynolds  number  based  on  the 

friction  velocity.  u„R/v0 
T  temperature 

7^  temperature  at  the  center 

T,  temperature  at  the  wall 

T,  characteristic  temperature,  q  c ,  put 

T.  dimensionless  temperature.  T  7, 

a.b  parameters  in  equations  ( 12H  14) 

cr  heat  capacity 

d  fixed-point  parameter  [5],  7  in  this 
work 

g°  bare  propagator  for  temperature 


q  heat  flux 

r,  wall  coordinate.  rum  \0 

u  mean  velocity  in  .v-direction 

u .  friction  velocity 

u,  dimensionless  velocity.  u  um 

r,  components  of  velocity. 

Greek  symbols 

a  inverse  total  Prandtl  number 

a0  inverse  molecular  Prandtl  number 

©  r-dependent  component  of 

temperature 

0.  dimensionless  ©,  Q,T, 

c  total  Prandtl  number,  v/k 

er0  molecular  Prandtl  number,  v0/k0 

e  expansion  parameter  in  RNG 

procedure 

i  turbulent  dissipation  rate 

k  total  diflusivity,  x0  +  K-lurb 

Ka  molecular  diffusivity 

/.  friction  coefficient,  8r,/pii,‘ 

/■o  expansion  parameter 

A,  integral  scale  of  turbulence 

p  fluid  density 

v  total  viscosity,  v0  +  vturb 

v0  molecular  viscosity 

v.  dimensionless  total  viscosity.  v  v0. 


outlined  in  the  Appendix.  The  details  of  the 
calculations  are  given  elsewhere  [5],  In  this  paper,  we 
are  interested  in  application  of  the  final  result  to  the 
problem  of  heat  transfer  in  turbulent  flow  in  a  pipe. 
This  will  be  done  in  the  next  section. 

The  most  distinguishing  characteristic  of  a 
turbulent  flow  is  approximate  universality  of  the 
properties  of  scales  much  smaller  than  any  integral 
scale  L  in  the  flow.  The  high  Reynolds  number 
turbulent  flow  is  characterized  by  three  different 
rangesofspatialscales.il  )  For  wave-numbers  L  >  n  L 
the  energy  spectrum  is  strongly  anisotropic  and  is  not 
characterized  in  any  universal  way  The  integral  scale 
reflects  both  geometry  of  the  flow  and  the  physico¬ 
chemical  processes  taktng  place  there.  (2)  At  much 
smaller  scales,  with  wave-numbers  satisfying  r.  L  <Z  k 
<  Lj  =  Re1  1 .  the  velocity  fluctuation  spectrum  is 
approximately  given  by  the  Kolmogorov  energy 
spectrum  f.tki  =  Cki:  r  with  the  Kolmogorov 

constant  C\  =  1  3-2.3.  (3)  In  the  dissipation  range 
(L  >  Ljl  the  energy  spectrum  decreases  exponentially 
with  k. 

Universality  of  the  small  scales  can  be  formulated  in 
the  language  of  theoretical  hydrodvn.imics  the  fluid 


described  by  the  Navier-Stokes  equation 
_  _£P  r:r, 

(  l  1  exj  r.v,  cx/  <  \j 
(V: 

-r1  =  0 
c.v, 


(4) 


subject  to  initial  and  boundary  conditions,  is 
characterized  at  the  small  scales  by  the  5  3- 
Kolmogorov  spectrum.  This  property  does  not 
depend  on  boundary  conditions  which  are  usually 
characterized  at  large  scales.  Boundary  conditions  can 
be  considered  from  the  viewpoint  of  small  scales  as  a 
source  of  energy  injected  into  the  large  scales  which 
subsequently  cascade  to  the  small  scales  Using  the 
analogy  with  equilibrium  statistical  mechanics  in 
which  the  results  are  independent  of  the  details  of  the 
interaction  of  the  system  with  a  heat  bath,  we  replace 
(4|  by  the  more  general  equation  in  and  add  i he  heat 
transfer  equation  |6| 


■  i , 

dr 


Heal  lr.in-.li-r  in  niihulent  fluids  I  Pipe  flow 


where  f  is  the  random  force  inotsei  chosen  to  generate 
the  velocity  Held  v  described  by  the  Kolmogorov 
spectrum  in  the  limn  of  large  wave-vectors  isrAll 
scales) 

It  has  been  shown  by  Yak  hot  [X]  that  the  Gaussian 
random  force  f  characterized  by  the  correlation 
function: 

< Jjk.vjytjik'.iy)}  ft  ck''*P1J(fc)i>(k+k')<>i(!j  +  i;/(  (7) 

with  P,j[k)  »  t\l-k,kj  k:.  generates  small-scale 
velocity  fluctuations  characterized  by  the  Kol¬ 
mogorov  spectrum.  The  parameter  t  in  (7) denotes  the 
dissipation  rate  of  the  turbulent  energy  per  unit  mass 
of  the  fluid  and  relates  the  force  f.  acting  on  small 
scales,  to  the  energy  input  taking  place  at  large  scales. 

This  fact  is  the  basis  for  using  the  random  force  <7| 
for  elimination  of  small  scales  in  the  construction  of 
either  turbulent  sub-grid  or  transport  models.  The 
renormalization-group  method  |RNG)  was  developed 
for  an  infinite,  homogeneous  medium  by  Forster  et  al. 
[6],  Martin  and  DeDominics  [7]  and  Yakhot  [8],  In 
these  works,  e  has  been  treated  as  a  given  parameter 
characterizing  the  rate  of  stirring.  In  finite  systems 


-  1  1  f  f 

e  =*  —  —  dr  1 6(.v.t)dx 

i*q  / cr,  ci-jY 
£<x.0  =  ~  • 

2  \cxj  cx,/ 


In  such  systems,  i  is  a  quantity  that  should  be 
determined  dynamically  from  the  equations  of  motion 
with  boundary  and  initial  conditions  applied.  The 
basic  ideas  of  the  renormalization  group  procedure  are 
given  in  the  Appendix. 

It  has  been  shown  by  Yakhot  and  Orszag  [9]  that 
the  Navier-Stokes  equations  for  the  mean  velocity 
field  r,  in  which  the  fluctuating  contributions  are 
removed  is: 

cV.  cv,  cp  c  cv, 

— +  ty— : — b-r~v  — .  (10) 

Ct  CXj  CX  ,  CXj  CXj 

Here  the  total  viscosity  v  takes  into  account  both 
molecular  and  turbulent  contributions  and  is  given  by 
the  following  relation  [9]: 


■M&-1 m)T 


where  the  ramp  function  H(x )  =  x  if  x  >  0  and 
Hix\  =  0  if  x  <  0  and  A,  is  the  inverse  integral  scale 
of  turbulence  [9],  The  parameter  = 

[2cflJ  +  2)]  =  0.333  since  if  =  7  for  this  problem  .  It  has 
been  shown  by  Yakhot  and  Orszag  [5]  that 
elimination  of  small  scales  from  the  equations  (4M6) 
of  a  passive  scalar  leads  to  the  following  relation 
between  thcinversetoial  Prandtl  number  x  =  a  ' 1  and 
the  total  viscosity  v: 


/  ’  -ii  ,  ,  x  -  />  . 1  l„ 

(,.. )L..J  ";i 

w  here 

_  a~  \ 

<i  -  h 

b  =  a  »-  I 

For  </  =  7.  relation  1 121  becomes 

/  x-  1.1793  Y>‘Y^:  I79-  V  " 

Vat,,  —  1.1793  )  V  x, ,  2. 1793  )  "  >  “4l 

The  result  (14|  expresses  the  inverse  total  Prandtl 
number  x  as  a  function  of  total  viscosity  i  and  is  the 
main  result  to  be  studied  in  this  paper. 

According  to  (II).  the  turbulent  viscosity  is  itself  a 
function  of  the  distance  from  the  wall  since  A,  must  be 
associated  with  the  distance  to  the  wall.  One  sees  that 
in  the  region  of  fully  developed  turbulence  where 
v0,  v  1.  the  total  Prandtl  number  a  =  x  ~ 1  =  0.8476. 
which  is  in  a  good  agreement  with  available 
experimental  data  a  =  0. 7-0.9  (see  Landau  and 
Lifshitz  [10]  and  Monin  and  Yaglom  [11]).  Close  to 
the  wall  where  v  %  v0,  one  finds  from  (14)  that 
a  ft  v0  k0.  Thus,  the  equation  of  motion  for  the  mean 
temperature  can  be  written  as: 

cT  ,cT  c  cT 

+  »—  k—  115) 

Ct  CX,  CX,  CX, 

where  k  -  avis  determined  from  (14).  The  dynamics  of 
diffusion  of  a  passive  scalar  is  governed  by  the  set  of 
equations  1 10),  ( 1 1).  ( 14).  ( 15). 

We  emphasize  that  these  results  do  not  include  any 
experimentally  adjustable  parameters. 


3.  HEAT  CONDUCTIVITY  IN  PIPE  FLOW 

Here  we  apply  the  results  presented  in  the  previous 
section  to  describe  the  process  of  heat  transfer  in 
turbulent  flow  through  a  pipe  of  radius  R.  The 
problem  can  be  formulated  in  terms  of  the  stationary 
Navier-Stokes  equation 


r  cr\  cr  j  cx 


and  the  heat  transfer  equation 

I r  /  cT\ 


I  r  /  cT 
-  —rx  — 
r  cr  \  cr 


where  v  and  k  are  total  viscosity  and  difTusivity, 
respectively.  The  parameters  v  and  x  include  b-'ih 
molecular  and  turbulent  contributions  The  total 
Prandtl  number  a  =  rv  is  determined  from 
relation  1 14). 

We  introduce  the  friction  velocity  u,  =  (t. 
wall  coordinate  r. ,  nondimension.il  velocity  «  . ,  and 


H«T  J0:l-* 
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nondimensional  total  viscosity  v. : 


r.  =r— ,  u.  =  uiu,,  v.  =  »7v0.  ( 16) 

vo 

The  equation  of  motion  now  has  the 
nondimensional  form: 

1  d  (  cu.\  2 

—  —  r.v*—  =  -—  (19) 

r.  cr,  \  cr,  )  R , 

where  Rm  «  u,R/v0  is  the  Reynolds  number  based  on 
the  friction  velocity. 

We  consider  heat  transfer  in  a  pipe  with  constant 
heat  flux  through  the  wall.  In  this  case  it  is  convenient 
to  introduce  a  new  variable  0  defined  as: 

T(x,r)*0(r)  +  Bx.  (20) 

Substituting  (20)  into  (17)  yields  an  equation  for  0(r): 


rrr  ^  cr  J 


t  . 

*  0  at  r  =  0 
cr 

(22) 

we  obtain: 

4a 

B  * - - — 

cppRe\0 

(23) 

where 

q  *  -ct,pK0)c©icr),,R 

(24) 

and  the  Reynolds  number 

2  r* 

Re  =  — -  u(r)r  dr. 

vo*  Jo 

(25) 

parameter  7*  is  defined  as  follows: 

T.  =  — —  .  (27) 

crl>u. 

Using  the  above  notation,  the  Nusselt  number  is 
given  by 

Su  =  C*o0Re  (28) 

where 

G  =  R;  1(2  f  ©.u.r.drA 


The  parameter  B  can  be  expressed  in  terms  of  the 
imposed  constant  heat  flux.  Integrating  (2 1 )  over  r  and 
using  the  fact  that 


Using  relations  (23M25I  the  heat  equation  (2 1 )  can  be 
written  in  the  nondimensional  form: 

1  c  /  .  \  2  u. 

—  —  — )*  — (-6) 
r*  rr„  \  cr  ^  J  Re 

where  0  .=0  1,  and  i  is  given  by  relation  ( 14).  The 


To  describe  heat  transfer  in  turbulence,  one  needs 
an  expression  for  the  coefficient  of  heat  conductivity 
which  takes  into  account  both  molecular  and 
turbulent  contributions  to  the  heat  transfer  process. 
The  theory  leading  to  relation  (14)  determines  the 
turbulent  diffusivity  in  terms  of  the  laminar  transport 
coefficients  and  the  turbulent  viscosity.  In  particular, 
it  describes  the  interaction  between  molecular  and 
turbulent  transport,  an  effect  of  much  significance  at 
low  Re  and  <j0.  Thus,  the  determination  of  turbulent 
heat  transfer  from  (14)  requires  reliable  data  on 
turbulent  viscosity.  Such  data  can  be  found  either  from 
theory  or  from  analysis  of  experimental  data  on 
velocity  profiles  in  pipe  flow. 

In  the  present  work  we  are  interested  exclusively  in 
demonstrating  the  power  of  the  ‘universal’  relation 
( 14)  provided  the  expression  for  turbulent  viscosity  is 
known.  Thus,  we  adopt  the  ad  hoc  model  [12]  for  the 
dimensionless  total  viscosity  v. : 

v.  =  l-t-0.41y.[l-exp(-is.  .4s)].  .4  =  26 

when  the  distance  to  the  wall  y,  <  50.  The  turbulent 
viscosity  for  y,  >  50  is  that  derived  from  the 
differential  k-l  mode!  of  Yakhot  and  Orszag  [5].  The 
model  viscosity  and  mean  velocity  profiles  obtained  by 
integrating  the  equation  of  motion  (19)  using  this 
viscosity  are  presented  in  Figs.  1  and  2.  The  friction 
coefficient  '/.  defined  by  r.  =  8  so 

/.  =  32|R.  Re)2  is  plotted  in  Fig.  3.  It  is  apparent  that 
the  agreement  with  experimental  data  is  very  good. 

The  equation  of  motion  ( 19i  and  heat  equation  1 26) 
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Flo  4  Dimensionless  temperature  profile  T.  -  T  1 
result*  of  calculation*.  experimental  Oatj  [  1 1] 


Fig.  2  Calculated  dimensionless  velocity  profile. 
11.  «  11  u..  Re  =-40.000. 


- 'Ju  =  c  018  Re 


Fto.  3.  Fnctton  coefficient  /.  lor  turbulent  pipe  flow:  x 
results  of  calculation  based  on  the  model  viscosity  from  Fig. 
1 : - Blasius  formula. 
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have  been  integrated  using  the  model  viscosity  from 
Fig.  1  and  a  from  the  relation  (14).  The  results  are 
presented  in  Figs.  4-10  for  various  Prandtl  (tr0)  and 
Reynolds  (R<?)  numbers. 

In  Fig.  4.  we  plot  the  calculated  and  measured 
temperature  profiles  for  air  flow  in  a  pipe.  As  we  can 
see  from  Fig.  4.  the  agreement  between  the 
experimental  data  and  the  results  of  calculations  for 
c0  =  0  7  is  very  good.  The  calculated  Nusselt  number 
for  air  flow  i<70  =  0  ~l  is  compared  in  Fig.  5  with  the 


Fig.  5.  Nusselt  number  Xu  as  a  function  of  Reynolds 
number  Re  =  u„£>  v0  for  the  air  flow  toc  =  0.7)  in  a  pipe  x 

results  of  calculations;  -  empirical  relation 

Xu  -  O.OI8RC0  V 


empirical  relation  widely  used  in  the  literature  [llj: 

Nu  =  0.018Re°  8.  (29) 


The  prediction  of  turbulent  heat  transfer  in  low- 
Prandti-number  flow  is  a  most  difficult  test  for  the 
model.  In  Fig.  6.  the  calculated  temperature  profiles  in 
liquid  mercury  (er0  =  0.02)  and  in  the  NaK  eutectic 
(a0  =  0.029)are  plotted  for  pipe  flow  at  Re  =  149.000 


x-  Rf.i49000  <r0  ■  0  7 
o-Bf., 49000  tr,  =0  029 
o- Re. 149000  CT0  *0  02 
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Fig  b  Dimensionless  temperature  1  T.  1  profiles  in  turbulent  flow  in  a  pipe  a!  Re  ■  140.000 
ia„  -  07|  NaK  cutest ic  M„  -  00291.  '  mercury  ia„  =  002), - from  ref  [I3J 
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Fig.  7.  Temperature  defect  (T.  -  T)  (T.  -  ^distribution  in 
turbulent  Dow  m  a  pipe:  x  results  of  calculation  for  NaK 
eutectic  (<r0  =  0.019); - experimental  data  of  ref.  [13]. 

It  is  apparent  that  some  fraction  of  the  temperature 
profile  (y.  =  100)  can  be  approximated  by  the 
logarithmic  law: 

T.  =  C  +4.85  log  y,  (30) 

where  the  constant  C  can  be  found  from  the  relation: 

C  =  3.1  +  3  In  <t0.  (31) 

A  relation  very  similar  to  (31)  has  been  obtained  from 
analysis  of  experimental  data  [11], 

It  should  be  mentioned  that  when  the  Prandtl 
number  is  small  (<r0  =  0.02)  the  logarithmic  part  of  the 
temperature  profile  appears  only  at  high  Reynolds 
numbers;  Re  >  10s.  For  Re  <  105,  no  part  of  the 
temperature  profile  can  be  approximated  by  the 
relation  (30). 

The  results  of  our  calculations  may  be  compared 
with  the  experimental  data  of  Buhr  et  al.  [13].  The 
temperature  profile  measured  in  the  NaK  eutectic  flow 
(<j0  «=  0.019)  at  Re  ^  4x  10*  is  compared  with  the 
present  results  in  Fig.  7.  In  Fig.  8.  we  plot  calculated 
and  experimental  data  for  the  temperature  profiles  in 
sexeral  low-Prandti-number  fluids  at  different 
Reynolds  numbers  in  the  range  3  x  10*  <  Re 
<  3.5  x  10s.  Again,  the  agreement  between  the  results 
of  calculations  and  experimental  data  is  very  good. 
The  Nusselt  number  .Yu  is  plotted  as  a  function  of  the 
Peclet  number  Pc  =  c0Rc  in  Fig.  9.  At  low  Pc.  the 
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Fig.  8.  Temperature  defect  (T»-T)  (T»  -  TC|  distribution  in 
turbulent  (low  in  a  pipe  as  a  function  of  Reynolds  number  Re: 
O  result  of  calculation; - experimental  data  of  ref.  [13]. 
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Fig.  9.  Nusselt  number  .Yu  as  a  function  of  Peclct  number 
Pc.  Su  s  6  74  when  Pc  5  0. 

results  of  numerical  calculations  gi\e  .Yu  %  6. "4.  This 
is  very  close  to  the  experimentally  obsened  [13.  14] 
limiting  Nussci.  number  Yu  '  6.S-'.0. 


I  h.  Mi  Niish-Ii  number  \u  .in  a  futieli.'n  of  Fr.imlll  number  r.  I  \petiriiem.ii  J.il.t  are  laxen  born 
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An.xhi  t  .Vni  of  hush  the  relation  l  I4i  and  of  the 
model  lor  turbulent  viscosity  jdopted  in  this  work  is 
the  prediction  of  heat  transfer  in  high-Prurultl-numbcr 
fluids  In  this  ease,  the  molecular  heat  diffusiv  its  is 
very  low  and  the  heat  transfer  process  is  determined 
entirely  by  the  turbulent  eddydiUusivits  The  results  of 
calculations  are  compared  with  experimental  data  in 
Fig  10  The  agreement  with  the  results  of 
measurements  [  1 1]  is  eery  good  across  a  wide  range  of 
Prandil  and  Reynolds  numbers.  I  <  a0  <  10”  and 
2  5  x  IQ4  <  Re  <  ;  x  10’. 

We  conclude  that  the  relation  (14)  can  be  used  for 
the  accurate  description  of  turbulent  heat  transfer 
throughout  the  entire  range  of  experimentally 
accessible  Prandtl  numbers,  which  vary  over  eight 
orders  of  magnitude,  i.e.  10' 1  <  tr0  <  lO". 
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Note  cutded  in  proof—  In  more  recent  work  on  the 
development  of  the  RNG  method  (see  ref.  [5]).  we  have  found 
that  the  proper  technique  is  to  evaluate  all  constants  at  the 
physical  dimension  </  =  3  to  lowest  order  in  an  expansion  in 
powers  of  c  rather  than  at  the  critical  dimension  d  »  7.  This 
modification  changes  the  turbulent  Prandtl  number  to  0.7 1 79 
and  changes  the  results  presented  here  by  several  percent. 
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APPENDIX 

In  this  Appendix  we  will  introduce  the  scale  elimination 
procedure  leading  to  renormalization  of  molecular  siscosity 
»■„  and  molecular  diffijsiviiy  n0.  Using  the  incomprcsstbiliiy 
condition  we  write  the  equations  of  motion  for  the  Fourier 


component*  nl  icJociiv  /.it  v/i.i/tj  passive  ssji.n  fit  <  o  is 
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1  *  d.j 
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2  J  i2x  i  '  ' 
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where  ./  is  the  dimensionality  of  the  space 

C'ltl  =  l -no- i,,t;i  1  i  A 3 1 

</°lt  I  m  I  —  KU  -  x,,i  ‘I  '  |.-\4| 

and  (he  random  force  f  i>  given  hy  ihe  correlation  function 
■  /,it  l((it  i  ■  *  i2."!|j  ' 1 2D.,t  /*  limit  -si  i  \  5 1 

The  projection  operator  P,n.  is  defined  as  P:.,iti  = 
t„  P„lk  I  t,  Plmlk  i.  Here  i  =  3  and  0  ,  /  /.. 

The  equations  of  motion  (All  and  i  A2i  are  defined  on  the 
domain  0  <  i  ^  A  The  RNG  procedure  consists  ol  mo 
steps.  First,  we  write  equations  in  terms  of  the  velocity  field 
decomposmon  or  two  components  i  Mi)  and  i  Mil  defined 
on  the  intervals  At"1  Sf  S  \  and  0  <  t  <  Ae  '. 
respectively  |/.0  «=  I ): 
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In  order  to  eliminate  modes  from  the  interval  Ae <  k  <  A. 
all  terms  i/lil.  T'lt)  should  be  removed  by  repeated 
substitution  of  1 A6)  for  v* .  T  *  back  into  i  A6 1  This  generates 
infinite  expansions  for  t*.T*  in  powers  of  /.0  in  which 
t >.T>  do  not  formally  appear.  Next,  averages  are  taken 
over  the  part  of  the  random  force  f  *  belonging  to  the  strip 
Ae*'  <  k  <  A.  This  procedure  formally  eliminates  the  modes 
Ae'1  <  k  <  A  from  the  problem. 

It  follows  from  (A6)  that,  after  removing  the  modes 
Ae'1  <  k  <  A.  the  equation  of  motion  for  r*.  T*  can  be 
written  up  to  second  order  in  J.0  as 
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The  parameter  d  *  7  at  the  fixed  point  Thus,  elimination  of 
small  scales  leads  to  renormalization  of  viscosity  and 
diffusivity.  The  second  step  of  the  procedure  consists  of 
iterating  the  scale-eliminatron  procedure.  This  leads  to  the 
results  given  in  Section  2. 


TRANSFERT  THERMIQUE  DANS  LES  ECOULEMENTS  TURBULENTS — 

1.  ECOULEMENT  DANS  UN  TUBE 

Resume — [.'expression  du  nombre  de  Prandtl  turbulent  obtenue  a  partir  d'une  procedure  de  groupe  de 
renormaltsalion  est  utilisee  pour  decrire  le  mecantsme  du  transfert  thermique  dans  1'ecoulement  turbulent 
dans  un  tube.  Les  resultats  sont  en  bon  accord  avec  des  donnees  expenmentales  dans  le  domatne  des 
nombres  de  Prandtl  10' :  <  o0  <  10*  accessibles  expenmentalemenl. 


wARMECBERGANG  IN  TURBULENTEN  FLUIDEN — I.  ROHRSTR0MUNG 

Zusammenfassung — Zur  Beschreibung  des  Warmeubergangs  bei  turbulenter  Rohrstromung  w-ird  der  Aus- 
druck  fur  die  turbulente  Prandtl-Zahl  verwendet.  den  man  aus  der  Renormalisations-Gruppen-Prozedur 
erhalt.  Die  Ergebmsse  stimmen  mil  experimentellen  Daten  tm  gesamten  Beretch  der  expenmentell  ver- 
fiigbaren  Prandtl-Zahlen.  10'*  <  <x0  <  10*.  gut  viberem. 


TEIUlOriEPEHOC  B  TyPEWIEHTHblX  5KH.3KOCTRX.  TEWEHHE  B  TPV’BE 

AjmoTaum — JLi*  onHcaHH*  TemonepeHoca  npit  Typoy.ieHTHOM  tcmchmh  ■  Tpvoe  Mcno.ibayeTca  supa- 
xeHHe  xia  TypoyneHTHoro  tHc.ta  flpaHam*.  no.ivHeHHoe  MeroaoM  peHopMaiHiauHOHHOH  rpynnu. 
Peav.TbTaTw  Haxoa*Tca  a  xopotneM  cooTeeTCTBHH  c  3KcnepHMeHTa.it.KUMM  p3hhumh  xis  Bcero  jHana- 
aona  3Ha4cHHit  HHC.ia  ripaHJT.i*,  10'J  <  £T0  <  lO6. 
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The  initiation  and  early  growth  of  apots  in  channel  and 
boundary  layer  flowa  ia  aimulated  uaing  a  three  dimensional 
apectral  code.  The  aimulated  apcta  show  significant 
agreement  with  available  experimental  data  for  such 
quantities  as  growth  rates  and  spreading  angles. 
Disturbances  are  introduced  into  the  center  and  edge  of  the 
developing  channel  spots  to  investigate  the  relative 
sensitivity  of  spots. 


INTRODUCTION 


1  . 

Emmons*  was  the  first  to  observe  turbulent  spots  in  a 
laminar  flow  undergoing  transition  to  a  turbulent  flow. 
Since  then  a  large  number  of  inves t igators  have  recognized 
the  importance  of  spots  in  the  study  of  both  transition  and 
turbulence.  Naturally  occurring  spots  are  initiated  by  flow 
disturbances  like  noise.  In  the  laboratory,  spots  may  be 
artificially  initiated  with  electric  sparks  or  by  injecting 
a  jet  of  fluid.  In  a  numerical  simulation  of  spots, 
controlled  disturbances  may  be  imposed  on  a  solution  of  the 
Navier-Stokes  equations. 

2 

Soon  after  Emnons '  discovery.  Elder  noted  that  spots 

tend  to  grow  independently  of  one  another,  even  when  they 

3  4 

overlap.  Caster  ’  studied  the  linear  growth  of  small 
amplitude  disturbances  into  a  wave  packet  using  both 
laboratory  experiments  and  theoretical  analysis.  His 
theoretical  predictions  have  been  confirmed  by  laboratory 
observations  so  long  as  nonlinear  effects  are  not  important. 
Wygnanski,  Sokolov,  and  Friedman^  conducted  an  experimental 
study  of  spots  in  a  boundary  layer.  Using  conditional 
sampling  techniques,  they  mapped  out  the  geometry  and  growth 
rates  of  a  spot  as  it  develops  in  a  boundary  layer.  Gad-el- 
Hak,  et  al.^  conducted  flow  visualization  experiments  on 
boundary  layer  spots  by  injecting  dye  upstream  of  the  spot 
initiation.  They  divided  the  spot  into  five  regionsCsee 
Figure  l).  Region  1  within  the  spot  overhangs  region  11, 


the  laminar  boundary  layer  below  the  head  of  the  spot. 
Region  III  appears  similar  to  a  turbulent  boundary  layer. 
In  regions  IV  and  V  the  flow  returns  to  a  "calm”  state.  The 
photograph  in  Figure  2  illustrates  the  characteristic 
arrowhead  shape  of  a  boundary  layer  spot  in  s t r eamwi s e - 
spanwise  projection.  This  photograph  was  obtained  by 
illuminating  dye  lines  with  a  sheet  of  light  very  close  to 
the  wall. 

The  first  detailed  research  directed  toward 

investigating  the  characteristics  of  spots  in  a  channel  was 
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conducted  by  Carlson,  et  al.  Using  mica  flakes  to 
visualize  the  flow  (Figure  3),  they  observed  that  a  channel 
spot  also  has  the  characteristic  arrowhead  shape.  They 
identified  (see  Figure  4)  several  features  present  in 
channel  spots.  The  spreading  half-angle  (l)  was  about  8 
degrees.  The  leading  edges  met  at  a  sharp  point  and  were 
preceded  by  oblique  waves(7).  The  center  of  the  spot  (4) 
contained  small  scale  turbulence.  Streaks(3)  trailed  from 
region  4. 

The  purpose  of  the  present  study  is  to  use  direct 
numerical  solution  of  the  Nav i e r - S t oke s  equation  to  identify 
details  of  the  internal  structure  of  spots,  as  well  as  to 
map  out  spot  dimensions  and  growth  rates.  Comparison  of  our 
results  for  growth  rates  of  the  large-scale  spot  dimensions 
with  those  seen  experimentally  verifies  that  the  essential 
growth  mechanisms  of  spots  is  captured  by  our  numerical 


One  previous  study  of  numerical  spots  should  be 
mentioned.  Leonard  used  discrete  vortex  methods  to 
simulate  numerically  the  early  growth  of  a  spot  ina 
boundary  layer.  As  with  the  present  computations,  the  spots 
computed  by  Leonard  are  typically  less  mature  than 
experimentally  observed  spots. 
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2.  OQMPUTAT I ONAL  GEOMETRIES  AND  NUMERICAL  METHODS 

The  computational  domain  that  we  use  to  simulate 
channel  flow  spots  is  as  follows.  In  our  simulations  of 
channel  flow  spots,  the  flow  is  represented  by  128x64x33 
Fourier  and  Chebyshev  modes  in  the  x  (streamwi se ) ,  y 
(spanwise),  and  z  (normal)  directions,  respectively  (see 
Figure  S).  The  flow  satisfies  periodic  boundary  conditions 
in  x  and  y  and  no-slip  (rigid)  boundary  conditions  at  the 
walls  (z-±l).  The  computational  box  is  nondimens  iona  1  i  zed 
by  the  channel  half-width;  in  the  runs  presented  below,  the 
physical  box  size  is  20x5x2.  With  128x64  resolution  in  x 
and  y,  the  resultant  node  spacing  (in  physical  space)  of  the 
spectral  collocation  points  is  4x»0.16  and  ^y-0.08. 

For  our  boundary  layer  spot  calculations,  the  flow  is 
represented  using  64  Fourier  modes  in  x  and  y,  with  Ai» 2  and 
jAy-1  (see  Figure  6).  In  the  z  direction,  the  33  collocation 
points  are  obtained  by  an  algebraic  mapping  of  the  interval 
[-1,1]  to  [O.oo]  with  half  the  collocation  points  located  in 
the  region  0<z<5.  The  computational  box  is 
nond  ime  ns  i  ona  1  i  ze  d  by  the  boundary  layer  thickness  f7-\/^xo  /U^ 
at  some  representative  x-location  xq .  The  periodic  boundary 
conditions  used  in  the  streamwise  direction  are  only 
approximate.  They  are  justified  because  the  increase  in 
boundary  layer  thickness  through  the  computational  domain  is 
only  6%  (see  also  Ba 1  a s ub r ama n i a n ,  et  al  ).  Wh  i  1  e  inflow- 
outflow  boundary  conditions  are,  in  principle,  more 
realistic  than  periodic  boundary  conditions,  they  are  more 


wasteful  of  spatial  resolution,  which  is  the  limiting  factor 
in  the  present  calculations. 

The  Nav i e r - S t oke s  equations  are  solved  in  rotational 

form, 

?.vxw  -  V(/7)  +  l/Re  V2(v) 
a  t 

V.  (v)-0 

where  w-Vxv  is  the  vorticity  and  J7-P+V  /2.  The  velocities 
are  normalized  with  respect  to  the  centerline  velocity  in 
the  channel  and  the  free  stream  velocity  in  the  boundary 
laye  r . 

The  spectral  method  of  Orszag  and  Patera10  is  used  in 
both  the  channel  and  boundary  layer  calculations.  For  the 
boundary  layer,  the  scheme  is  modified  by  mapping  the 
Chebyshev  collocation  points  of  the  channel  to  the  desired 
locations  in  the  boundary  layer.  A  mapping  function 

z*  -  f (z) 

is  chosen.  When  taking  derivatives  in  the  z-direction 
(e.g.,  in  calculating  the  vorticity)  the  Chebyshev 
differentiation  in  z*  is  followed  by  multiplication  by 


The  boundary  condition  at  infinity  (v^  -  1)  is  implemented 
by  recalculating  the  (0,0)  Fourier  mode  (the  mean  flow  in  z 
and  y)  in  the  viscous  step.  Symmetry  is  not  imposed,  but 
the  spots  develop  symmetrically  when  symmetric  initial 
conditions  are  used. 

The  disturbance  is  initiated  by  applying  a  body  force 
to  a  packet  of  fluid,  producing  a  small  jet  normal  to  and 
away  from  the  wall.  The  form  of  the  disturbance  is  Gaussian 
in  z  and  y  and  continuous  in  time. 


F 


G(t) 


t[-r2/2c2) 


where  G(t)  is  a  ramp  function.  The  size  of  the  jets  are 
indicated  in  Table  I. 


table  I 

Channel  B.L. 


o  0.16 

Location  0.1  <  z  <  0.2 
Peak  normal 

Velocity  0.09 


0.7 

0.05  <  z  <  1.5 


0.035 


We  impose  the  following  boundary  conditions  on  the  flow 
through  the  channel:  the  velocity  at  the  walls  is  zero,  and 
the  flow  is  periodic  at  the  inflow/outflow  and  cross- stream 


boundaries.  The  Reynolds  number  for  the  channel  runs  is 
6000  based  on  the  channel  half-width.  The  Reynolds  number 
for  the  boundary  layer  simulations  is  1000  based  on  the 
boundary  layer  thickness  corresponding  to  r?« 1.0. 

3.  SPOTS  IN  CHANNEL  FLCWS 

In  Figure  7,  we  plot  contours  of  the  maximum  (in  y)  of 

the  absolute  value  of  the  normal  velocity,  Max  Iv  I  for  a 

y  z 

channel  spot  at  Re-6000.  The  contour  plots  we  present  for 
channel  spots  encompass  the  entire  20x5x2  computational 
domain;  their  dimensions  are  not  to  scale.  Except  where 
noted,  the  contours  are  at  1%  intervals  of  Max  Iv  I,  where  q 
is  the  coordinate  normal  to  the  plotting  plane.  With  this 
projection  of  the  spot  onto  a  plane  we  view  the  data  from 
the  experimentalist’s  perspective  (with  the  line  of  sight 
extending  all  the  way  through  the  channel).  At  time  t-1 , 
the  initial  disturbance  has  convected  downstream  and  has 
become  slightly  distorted.  The  initial  peak  velocity  of 
0.09  has  decreased  to  0.038  due  to  viscous  diffusion.  By  a 
time  of  t-3,  the  velocity  has  increased  to  0.05  due  to 
instability  in  the  flow.  The  disturbance  is  elongated  in  x 
as  well  as  convected  downstream.  In  Figure  8  we  see  the 
disturbance  develop  most  of  the  features  characteristic  of  a 
spot.  The  front  of  the  disturbance  moves  away  from  the 
wall.  The  disturbance  grows  in  all  directions  and  the 
"arrowhead"  shape  becomes  apparent.  The  peak  normal 

velocity  increases  from  670  at  t  —  1 2  to  99b  at  t-18. 


In  Figure  9  we  show  the  development  of  the  boundaries 
with  an  isometric  view.  Enclosed  within  the  surface  is 
fluid  whose  x  velocity  differs  from  the  Poiseuille  profile 
by  more  than  2%. 

The  results  plotted  in  Figure  10  are  Max  Iv  I  and 

y  ® 

Max^Iv^l.  At  t-30,  the  spot  has  fully  extended  through  the 
channel  with  a  peak  normal  velocity  of  13%.  The  initial 
disturbance  on  the  bottom  wall  has  induced  a  new  disturbance 
at  the  top  wall.  This  second,  smaller  spot  has  a  peak 
velocity  that  occurs  at  a  distance  of  approximately  0.25 
(1/8  channel  width)  away  from  the  top  wall.  By  t-30,  the 
two  spots  have  joined  to  produce  a  disturbance  that  fills 
the  span  of  the  channel. 

In  Figure  11  we  show  the  distortion  of  the  Poiseuille 
profile  at  the  spot  center.  The  velocity  at  the  edge  is 
essentially  unchanged  from  that  of  the  original  Poiseuille 
flow,  while  at  the  spot  center  there  is  a  velocity  defect  of 
0. 1-0.2.  At  the  bottom  wall,  the  shear  has  increased  by  a 
factor  of  3 . 

In  Table  II  and  Figure  12,  we  show  how  the  spot 

geometry  changes  in  time.  Although  there  are  significant 

differences  between  conditions  generating  our  numerical  spot 

and  those  generating  the  spots  studied  experimentally,  a 

comparison  of  numerical  and  laboratory  features  is 

8 

instructive.  Carlson  et  al  generated  spots  in  a  laboratory 
channel  flow  at  Re-1000,  while  we  used  Re-6000  in  our 
calculations.  Most  of  the  experimental  data  were  taken  more 


than  50-100  channel  widths  downstream  of  the  initial 
disturbance,  while  we  have  been  able  to  follow  the  spot  for 
only  10  channel  widths.  Further  development  of  the  channel 
spot  would  require  a  larger  computational  domain.  The 
growth  rate  of  the  width  and  length  of  the  numerical  spot 
becomes  constant  at  t«15  and  remains  so  until  the  spot  fills 
the  domain  at  t-32.  This  steady  growth  rate  is  slightly 
higher  than  that  observed  experimentally  in  both  the  lateral 
and  longitudinal  directions.  This  discrepancy  can  be  due  to 
the  difference  in  Reynolds  numbers  or  to  the  lack,  of 
maturity  of  our  computed  spots  compared  to  those  studied 
experimentally.  We  have  not  observed  in  our  data  any 
significant  evidence  of  the  leading  To  1 lmi en-Schl i c t ing 
waves  that  were  observed  experimentally.  Again,  we  believe 
that  the  absence  of  these  waves  is  due  to  the  lack  of 
maturity  of  our  computed  spots. 

Table  II  Channel  Spots 

Experimental  I  Computational 


Velocity  of  Front  0.6  I  0.85 

Rear  0.34  I  0.25 

I 

8°  I  10° 


Spreading  Half-Angle 


A  further  numerical  calculation  was  done  to  compare  the 
stability  characteristics  of  the  spot  at  its  edge  and 
center.  The  velocity  field  of  a  spot  at  t-20  is  used  as  the 
initial  condition  for  three  runs.  The  first  run  consists  of 
the  restarting  the  original  spot  calculation  and  allowing 
the  spot  to  continue  development  undisturbed  to  a  time  of 
24.  For  the  second  run,  a  disturbance  is  applied  at  t-20  to 
the  original  spot  at  its  center.  This  disturbance  is  of  the 
same  spatial  and  temporal  extent  as  the  original  disturbance 
that  initiated  the  spot,  but  the  magnitude  is  1/10  that  of 
the  original.  The  difference  between  the  resulting  velocity 
fields,  e(x,t)  -  I  vzl  *VZ2  *•  4  measure  of  the  effect  of 
the  disturbance.  By  t-24,  e(x,t)  has  exceeded  1%  in  the 
central  2/3  of  the  spot  (Figure  13).  The  third  run  is 
identical  to  the  second,  but  with  the  disturbance  applied  at 
the  spot  edge,  rather  than  at  the  center.  At  t-24  the 
disturbance  had  propagated  through  mos  t  of  the  spot  (see 
Figure  14),  and  had  a  peak,  magnitude  of  about  4%,  as  opposed 
to  the  1.5%  peak  from  the  second  run. 

From  these  results,  we  conclude  that  channel  spots  are 
less  stable  at  their  edges  than  at  their  centers.  This 
observation  suggests  that  spots  grow  by  destabilization  of 
neighboring  fluid,  rather  than  simply  engulfing  laminar 
fluid. 


4.  SPOTS  IN  BOUNDARY  LAYER  FLOWS 


2 


In  Figures  15  through  17  we  show  the  development  of  a 

boundary  layer  spot  at  Re-1000  up  to  t-90.  The  contour  plots 

we  present  for  boundary  layer  spots  encompass  the  entire 

128x64  computational  domain  in  x  and  y  and  are  truncated  at 

z-22.  Again,  except  where  noted,  the  contours  are  at  1% 

intervals  of  Max  Iv  1,  where  q  is  the  direction  normal  to 

q  z 

the  plotting  plane.  Figure  15  shows  the  streamwise  and 
spanwise  development  of  the  spot  from  the  initial 
disturbance.  At  t-90,  the  spot  has  begun  to  develop  the 
characteristic  arrowhead  shape,  which  is  more  apparent  in 
the  second  (29b)  velocity  contour.  Figure  16  shows  the 
development  of  the  triangular  shape  and  the  overhang  in  the 
spanwise  direction.  Figure  17  shows  the  overhang  develop  in 
the  leading  edge. 


Table  III  Boundary  Layer  Spots 


Experimental  I  Computational 


Velocity  of  Front  0.9 

1 

0.85 

Rear  0.5 

1 

1 

0.3 

Spreading  Half-Angle  10° 

1 

1 

12° 

The  growth  and  development 

of  the  spot 

i n  a  bounds ry 

layer  is  compared  with  the 

experimental 

findings  of 

Wygnanski,  et  al^  in  Table  III.  The  growth  rate  of  the  spot 


in  the  streamwise  and  spanwise  directions  is  in  relatively 
close  agreement  with  the  experimental  data.  This  suggests 
that  the  growth  mechanisms  in  a  boundary  layer  spot  have 
been  accurately  captured  in  this  simulation. 


Figure 

18  shows  cross 

sections  of  the 

spot 

a  t 

t-90. 

Here  we  plot 

contours  of 

the 

local 

values 

of 

vz 

at  y- 

y  .  -0.5, 

J center 

2.5,  and  4.5, 

in 

Figures 

18a, 

18b, 

and 

18c  , 

respect ive ly 

.  I n t e rva Is  are 

a  t 

1%  and 

dashed 

contour 

lines 

represent  negative  z  velocities.  The  velocities  are  highest 
in  the  plane  closest  to  the  center  of  the  spot  (see  Figure 
18a).  Away  from  the  spot  centerline  the  velocities  and  the 
spot  height  decrease.  The  front  of  the  spot  has  an  overhang 
of  a  distance  of  10-20  in  x,  as  has  been  observed 
experimentally.  The  flow  is  dominated  by  eddies  with  length 
scales  of  approximately  10  in  x  and  5  in  y.  These  length 
scales  differ  from  those  of  unstable  modes  of  the  Orr- 
Sonmerfeld  equations,  which  predicts  linear  instability  for 
much  longer  wavelengths,  30<)vx<85. 

In  order  to  explore  the  later  time  evolution  of 
boundary  layer  spots,  it  will  be  necessary  to  use  higher 
resolution  simulations,  which  we  hope  to  perform  in  the 
future  . 

CONCLUSIONS 

It  has  been  sh own  that  spots  can  be  generated  by 
numerical  solution  of  the  Navier-Stok.es  equations.  The  fact 
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that  our  results  for  the  growth  rates  of  the  large-scale 
•  pot  d  imens  i  ons  are  relatively  close  to  those  seen 
experimentally  suggests  that  the  essential  growth  mechani sms 
of  spots  have  been  captured  by  our  numerical  experiments. 
These  simulated  spots  are  less  mature  than  typical 
experimental  spots,  but  their  behavior  appears  to 
approximate  that  in  a  fully  developed  spot. 

The  spots  generated  were  not  dominated  by  two 
dimensional  To  1 lmi  en- Sch 1 i c t ing  wave s  .  This  suggests  that 
the  growth  in  spots  is  not  linear  growth  of  two  dimensional 
To  1 lmi  en-Schl ict ing  waves.  Moreover,  the  perturbation 
velocities  seen  were  about  0.1;  perturbations  this  large 
would  make  the  results  of  linear  theory  inapplicable  and 
suggest  domination  of  nonlinear  effects.  This  does  not  rule 
out  the  importance  of  To  1 lmi  en- Schl i c t ing  waves  in  the 
amplification  of  small  disturbances  which  may  develop  into 
spots  or  as  a  driving  mechanism  for  some  secondary 


instability  in  spots. 
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FIG.  1  Schematic  of  an  experimental  boundary  layer  spot  cut 
through  the  center  ( f  rom  Gad-e  1  -Hak  et  al.^). 

FIG.  2  Visualization  of  an  experimental  boundary  layer  spot 
using  fluorescent  dye  and  a  sheet  of  laser  light  at  the 
wall;  Rex-5xl0^  (f rom  Gad-e 1 -Hak  et  al.^). 

FIG.  3  Visualization  of  an  experimental  channel  spot  using 

g 

mica  platelets  (from  Carlson  et  al.  ). 

FIG.  4  Channel  spot  schematic:  (l)  spreading  half  angle;  (2) 
trailing  streaks;  (3)  region  of  small-scale  turbulence  (4) 
oblique  Tol lmi  en-Schl i c t ing  waves  (from  Carlson  et  al.  ). 
FIG.  5  Channel  geometry  and  nomenclature.  Channel  is  20x5x2 
in  the  x,y.  and  z  directions,  with  128x64  Fourier  modes  in  x 
and  y  and  33  Chebyshev  modes  in  z. 

FIG.  6  Boundary  layer  geometry  and  nomenclature.  Boundary 
layer  computational  domain  is  128x64  in  the  x  and  y 
directions,  with  64x64  Fourier  modes  in  x  and  y  and  33 
Chebyshev  modes  mapped  in  the  normal(z)  direction. 

FIG.  7  Early-time  evolution  of  channel  spot.  Max  Iv  I 

y  * 

contours  are  plotted  at  1%  intervals. 

FIG.  8  Channel  spot  at  intermediate  times.  Maxzlvzl 
contours  in  a)  and  b);  Max  Iv  I  contours  in  c)  and  d). 

y  * 

FIG.  9  Surfaces  of  2%  x-velocity  perturbations  in  developing 
channe  1  spot  . 

FIG.  10  Channel  spot  at  t<*30.  Maxzlvzl  contours  in  a); 
Max  Iv  I  contours  in  b). 

y  z 

FIG.  11  Mean  velocity  profiles  at  center  (solid)  and  edge 
(broken)  of  spot . 
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FIG.  12  Location  in  x  of  front,  center,  and  rear  of  channel 
•pot  vs.  time,  where  spot  is  defined  as  region  where 
lvzl^2%.  For  t  larger  than  30,  the  spot  length  reaches  the 
periodicity  length  of  the  computational  domain,  so  the  spot 
ceases  to  grow  in  the  streamwise  direction. 

FIG.  13  Perturbation  velocity,  e(x,t),  contours  at  t-22  and 

t-24  for  channel  spot  perturbed  at  its  center  at  t-20. 

FIG.  14  Perturbation  velocity,  <(x,t),  contours  at  t-22  and 

t-24  for  channel  spot  perturbed  at  its  edge  at  t-20. 

FIG.  15  Development  of  boundary  layer  spot.  Max  Iv  I 

z  z 

contours  are  plotted  at  1%  intervals.  a)t-30;  b)t-60;  c)t-90 

FIG.  16  Development  of  boundary  layer  spot.  Max  Iv  I 

y  z 

contours  are  plotted  at  1%  intervals.  a)t-30;  b)t-60;  c)t-90 

FIG.  17  Development  of  boundary  layer  spot.  Max  Iv  I 

z  z 

contours  are  plotted  at  1%  intervals.  a)t-30;  b)t-60;  c)t-90 

FIG.  18  Slices  of  spot  at  t-90.  Contours  of  v  at  y-32,  30, 

z 

and  28.  The  plane  of  symmetry  of  the  spot  is  at  y-32. 5. 
Dotted  lines  represent  negative  v^ . 
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Figure 


Early-Time  Spot  Evolution  at  R 
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Figure  7 


MEAN  VELOCITY  PROFILES 


AT  CENTER  AND  EDGE  OF  SPOT 
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